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Abstract 

Extending the Swendsen-Wang cluster algorithm to include both bulk (H) 
and surface fields (Hi ) in L x L x D Ising films of thickness D and two free LxL 
surfaces, a Monte Carlo study of the capillary condensation critical point of 
the model is presented. Applying a finite-size scaling analysis where the lateral 
linear dimension L is varied over a wide range, the critical temperature T C (D) 
and the associated critical field H C (D) are estimated for 4 < D < 32 lattice 
spacings, for a choice of the surface field H\ small enough that the dependence 
of H C (D) on Hi is still linear. It is shown that the results are consistent 
with the power laws predicted by Fisher and Nakanishi [M.E. Fisher and H. 
Nakanishi, J. Chem. Phys. 75, 5857 (1981)], namely T c (oo)-T c (D) oc D~ x l v , 
H C (D) oc Z)-(A-Ai)/^ where 

v is the bulk correlation length exponent of 
the three-dimensional Ising model, and A, Ai are the corresponding "gap 
exponents" associated with bulk and surface fields, respectively. As expected, 
the order parameter of the thin film near its critical point exhibits critical 
behavior compatible with the universality class of the two-dimensional Ising 
model. 



I. INTRODUCTION 



The application prospects of nanoscale technology have created a fresh interest in the 
behavior of both simple fluids and complex fluids confined in pores or in a thin film geometry 
in layers confined by parallel wallJ~§. However, a prerequisite for the clarification of pattern 
formationill and dynamics@i is a good understanding of the interplay between bulk and sur- 
face effects on thermodynamics and phase behavior in this finite-size geometryiHl: although 
theoretical aspects of phase transitions and critical phenomena in confined geometry have 
been considered since a long time&^l, this is still a topic of active current researchiHIl 
even for one of the most well-known phenomena, namely "capillary condensation" 0H. By 
"capillary condensation" one means the finding, already discovered in the 19th centuryS, 
that in a capillary the condensation of a gas occurs already at a lower pressure p than the 
coexistence pressure p coe ^ necessary to induce condensation in the bulk. Qualitatively, this 
shift of the transition can be attributed to the interaction of the fluid molecules with the 
attractive walls of the capillary. Although confinement effects on fluids and their phase 
transitions have been studied experimentally since a long time as welllHil, a quantitative 
characterization of the shift of the capillary condensation critical point remains a challenge. 
While for temperatures T below the critical temperature T C (D) of the thin film of thickness 
D the chemical potential at the condensation transition p, c {D) is shifted relative to its bulk 
value simply as ^ C {D) — /U c (oo) oc D^ 1 ("Kelvin equation" )0, for large enough D, Fisher 
and NakanishiEl predicted a completely different power law for the corresponding shift at 
the critical temperature itself, namely 

/x c (D) - /x c (oo) oc IT ( A - A ^ T = T C {D) (1) 

for weak surface forces. In Eq. ([I]), critical exponents of the three-dimensional Ising model 
universality class (that encompasses criticality of gas-fluid critical points or the related un- 
mixing transitions in binary mixtures, etc.) enter, namely the correlation length exponent^ 
v ~ 0.63 and the "gap exponent" A = 7 + (3 ~ 1.56 and the corresponding exponent for a 
free surface0'0~H Ai ~ 0.46 — 0.48. Also for the shift of T c a similar power law holds, 



T c (oo) - T C (D) oc D- 1 /", (2) 

which is the same relation as is familiar from standard finite-size scalingSllIHil for the shift 
of T c in films with "neutral walls" {i.e., no surface field preferring one of the phases coexisting 
for T < T C (D) at fJ> c {D) = [i c (oo) act} or in films with periodic boundary conditions where 
surface effects are a priori absent. While for the latter systems Eq. (Q) has been studied 
by various methods0'0&~0, Eq. ([!]) to our knowledge has not yet been tested by Monte 
Carlo simulations. In previous workHSil, tests of the Kelvin equation and corrections to 
it@ have been carried out and a capillary condensation critical point was located for a thin 
Ising film@ but for a single value of D only. 

In the present paper, we fill in this gap by presenting a Monte Carlo study of the critical 
behavior of capillary condensation in thin Ising films for a range of thicknesses. Invoking 
the universality principle!!, one can argue that nearest neighbor Ising lattices with short 
range surface fields should yield the same power law, Eq. (|l|), as more realistic models and 
real fluids in slit-like capillaries do. Unlike the situation in real fluids, packing effects at 
the surfaces and a dependence of the density in the middle of the film on its thickness are, 
however, absent and one focuses on the universal critical behavior. One can argue that the 
order parameter correlations in the directions parallel to the wall should all scale with the 
critical exponents of the universality class of the two-dimensional Ising modeS, 

MocF 2 , xoc|t|" 72 , £|| oc \t\-»>, #8 = 1/8, 72 = ^, v 2 = 1, 

t=l- T/T C (D) -> 0, all D < oo. (3) 

Of course, one expects that for large D the asymptotic critical region where this 
two-dimensional critical behavior holds is very narrow, due to a crossover to the three- 
dimensional critical behavior as D — > oo, and a quantitative understanding of this 
crossover0'il~il is a challenging aspect of this problem, too. 

In Sec. II, we shall hence briefly define the model that is studied and the quantities 
that will be analyzed and comment on the simulation methods. Sec. Ill briefly reviews the 



scaling predictions, including the finite-size scaling results for the case where both D and 
the lateral linear dimension L are finite. Sec. IV then presents our results on T C (D) and the 
critical fields H C (D), on which our tests of Eqs. ([[]), (0) are based. Sec. V discusses those 
aspects of our results which are pertinent to a test of two-dimensional criticality, Eq. ([3|), 
while Sec. VI summarizes our conclusions. 



II. MODEL AND SIMULATION TECHNIQUE 

Invoking the well-known isomorphism between the lattice gas model of fluids and the 
Ising model of magnetism (see e.g. Ref. ^ for details), we study the Ising model on the 
simple cubic lattice in the presence of a bulk field H and a surface field Hi, 

H = -jJ2S i S J -Hj2S l -H 1 Y, s » S i = ±l, (4) 

(ij) i it surfaces 

where the exchange interaction J is only present between nearest neighbors on the lattice. 
Note that phase coexistence in the bulk (phases with positive and negative magnetization 
correspond to gas and liquid phases of the fluid, respectively) corresponds to H = 0, see 
Fig. [I]. In the thin film, one trivially obtains the result that for zero temperature phase coex- 
istence occurs fori! H cocx (D, T = 0) = —2Hi/D, but for T > the variation of H cocx (D, T) 
is nontrivial. While previous workU was mostly interested in the behavior of H cocx (D,T) 
near the temperature T w (Hi) of the wetting transitioni'lHH, we consider here scaled sur- 
face fields HiD Al / u small enough such that we stay in the nonwet regime of the surface 
phase diagram of the semi-infinite systemEl throughout, although we consider the vicinity 
of T c °° = T C (D = oo). Measuring all lengths in units of the lattice spacing, we consider film 
thicknesses D = 4, 8, 12, 16, 24, 28 and 32 for an L x L x D geometry, varying L over an as 
wide range as is practical, from the point of view of available computer ressources. In the 
x, y-directions parallel to the thin film, we apply periodic boundary conditions as usual0@. 

In order to be able to find T C (D) and H C (D) = H coex (D,T = T C (D)) reliably, we have 
to use aspect ratios L/D ^> 1. Although the choices of film thickness as quoted above are 
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not extremely large, it is clear that use of rather large linear dimensions L is mandatory 
for obtaining reliable results. If we would use the Metropolis algorithmIl~ll, as done in 
Refs. |T3|J2"8| , or the heatbath algorithm^, "critical slowing down" 00 would be a serious 
problem: i.e., the "time" r over which subsequently generated system configurations are 
correlated varies as0 

r oc L* with z « 2.16 (d = 2) or z « 2.03 (d = 3), (5) 

the prefactor in this power law being of order unity if Monte Carlo time is measured in units 
of attempted Monte Carlo steps (MCS) per spin. Since we wish to use linear dimensions of 
the order of L 10 2 , relaxation times of the order of 10 4 MCS easily result. Given the fact 
that quantities like the specific heat C v and the susceptibility x, recorded from fluctuations 
of energy and magnetization, 

C v = ((h 2 ) - (H) 2 ) I (L 2 Dk B T 2 ), x = j / (E 5 <) ) - (E 5 *) } /( L2Dk BT), (6) 



are non-self-averagingEH^Ea, one needs n ^> 1 statistically independent observations (i.e., 
separated by time intervals At > r) to obtain C v and x with small enough error (the 



relative error of these quantities isS y 2/n, irrespective of L and D). For this reason, it is 
clear that the use of cluster algorithms which reduce critical slowing down0IHll is highly 
desirable. However, for the present problem where both a bulk magnetic field and a surface 
magnetic field of competing sign are present {Eq. (^)} application of cluster algorithms is 
nontrivial. It turns out that an extension of the "ghost spin algorithm" to the present 
problem is rather straightforwardly possible@@. The coupling of spins to a magnetic field 
is thereby treated as if it were an additional infinite-range exchange coupling to an extra 
spin Sg = ±1- This coupling has the strength h = \H\ for spins in the interior of the film 
and h — \Hi + H\ for spins in the surface layers. In addition to putting bonds in clusters of 
spins (inside a cluster all spins are connected by exchange interactions and have the same 
sign) with probability^^ pb = 1 — exp(— 2J/k B T) one also puts bonds between the spins 
in clusters and the ghost spin pc = 1 — exp(— 2hjk B T\ if the orientation of the spins in the 



cluster is the same as that of the ghost spin {which is Sq =sign (H) for interior spins and 
Sq =sign (Hi + H) for spins in the surface planes, respectively}. 

While this extension of the cluster algorithm to the case of nonzero bulk and surface 
fields is formally exact, discussion of its efficiency is a rather delicate problem: in fact, if 
h/k^T is of order unity, also pa is of order unity and the infinite-range character of this 
coupling then implies that huge clusters containing a large fraction of the entire simulation 
volume would be created most of the time! It is clear that under such circumstances the 
algorithm would be very inefficient; as in the case of zero field it is necessary for a good 
performance of a cluster algorithm that typically large clusters are created but a single large 
cluster must contain only a negligible fraction of the total volume in the thermodynamic 
limit. As a consequence, one needs h/k^T <C 1, and since k-^T is in the range of 3.5 — 4.5 
we have thus chosen to work with a single value of the surface field, namely Hi = —0.015 J. 
Even for this small value - note that the corresponding value of H is typically one or two 
orders of magnitude smaller, see below - the performance of the algorithm has significantly 
deteriorated, in comparison with the case without any magnetic fields. This fact can be 
clearly demonstrated by a binning analysis^lS of the magnetization m in the system: the 

N 

N (dynamically correlated) subsequent observations m u = (1/N)J2S^ are grouped into 

i=i 

n = N/Nt blocks of length Nb, from which block averages m^ of the corresponding Nb 
observations {m u } belonging to the block with index \i are formed. Then 



Am = [n(n - 1)]~ 1/2 , 



\ JI(m M -m) 2 (7) 
\ i=i 

(where m = rT l J2 m i) is studied as a function of block length Nb (Fig. [|): when Am 

8=1 

is independent of Nb, the subsequent m^ are statistically independent, and Am is a good 
estimation of the statistical error; otherwise one sees a systematic increase of Am with Nb 
and the value of Nb needed to reach a saturation value yields an estimate for the correlation 
time. For the Metropolis algorithm and the chosen system size (L=128), even for Nb = 4000 
one is far from saturation, and hence it is clear that this algorithm would be very impractical 
for the present problem. For the cluster algorithm and Hi = 0, on the other hand, Am 



vs. Nb is essentially constant, Am pa 0.0004, the correlation time being of order unity, as 
expectedH'il'0lll. However, this is not so for the cluster algorithm in the case Hi = —0.015: 
Am saturates at a plateau of about Am pa 0.007, i.e. the error is almost a factor 20 larger, 
and the correlation time is of the order of r m pa 280 Monte Carlo steps in the example shown 
in Fig. 0(b). Thus, while the gain of the cluster algorithm in the zero field case compared to 
the Metropolis algorithm is very significant, in our problem it is only rather modest! This 
- somewhat unexpected - dramatic decrease of the efficiency of the cluster algorithm with 
increasing strength of the surface (and bulk) fields has prevented us both from studying 
systems larger than L = 128 and from studying the dependence on Hi systematically. Runs 
of length up to 1.2 million Monte Carlo steps (MCS) were performed. 

As is well knownllElzl, cluster algorithms at critical points of Ising systems are rather 
sensitive to correlations among the pseudorandom numbers produced by the random number 
generator. In the present work, we used thus an improved version of the standard "R250" 
generator^, where two versions {one based on the pair of integers (250, 103) the other with 
the pair (521, 168)} are combined with the logical XOR operation. 

In order to make best use of our simulation data, we apply standard multihistogram 
interpolation techniques00ElIll. Note that we needed a three-dimensional histogram 
P(E,m,mi), E being the exchange energy, mi the magnetization in the surface plane, 
in order to allow reweightings in the full parameter space of independent control variables 
(T, H, Hi) and hence the storage requirements for P are nontrivial. However, noting that 
all measurements of E, m, mi can be represented by integers, each integer needing 4 Byte, 
we can store the time series of 10 6 observations with a storage of 12 MByte, irrespective of 
the choices of L and D. 

The multihistogram reweighting with respect to the bulk field H is crucial in order to 
be able to find the field if coex (T), along which for T < T C (D) two-phase coexistence occurs, 
applying the "equal weight rule"S0@§ In the space of variables (£7, m, mi), the two 
phases show up as separate peaks of P(E,m,mi) {or P(m, m^, respectively, see Fig. |3](a), 
when one studies an isotherm one can integrate out E, of course}, which have precisely 
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the same weight at H = H coex (T) while for H ^ if coex (T) (but not too far away from it) 
the two peaks can still be identified but have different weight. With the multihistogram 
reweighting, a small number of simulation points suffices to generate the curve H coex (T) 
{and its extrapolation into the regime T > T C (D)} with reasonable precision, see Fig. 0(b). 
Since near T C (D) the free energy differences between the two phases are very small also off 
coexistence, the statistical error in the estimation of H coex (T) is not negligible, and also 
systematic errors, since L/D is not large enough, need to be considered. The latter problem 
also affects the estimation of T C (D), as will be discussed in Sec. IV. 

III. SCALING PREDICTIONS 

For completeness, we first summarize the pertinent predictions of the scaling theory 
for thin Ising films near the critical point&§S, assuming the lateral linear dimension L 
infinite, and consider the extension!! to finite L in the following. The singular part of the 
free energy per spin is assumed to scale as follows 

/sing (D,T,H,Hi) « \t\ 2 ~ a f ± (D\t\ v , H\t\~ A , tfi|C Al ) , (8) 

where a is the exponent of the specific heat of the three-dimensional Ising model, t = 
(T — T c (oo)) /T c (oo) , f± is a "scaling function" (with two different branches, referring to 
the sign of t), and the other exponents have already been defined in Sec. I. 
Now it is convenient to introduce the scaling variables 

x = D\t\" , w = H x D^l v , (9) 

and then Eq. (H) can also be written as, eliminating \t\ from the arguments of f±, 

/sing (D, T, H, ITO » \tr«f ± (x, , . (10) 

Since the critical point of the thin film is shifted relative to the bulk critical point T c (oo), 
it must correspond to a singular behavior of the scaling function f±. At fixed Hi and fixed D 
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this means the scaling function f± (x,y,yi = w/x Al ^^j has a singularity at a point x c (w) , 
y c (tw). Therefore the shifts AT C (£>) , AH C (D) follow as§ 

AT C = T c (D, H 1 ) - T c (oo) = -B T D~ l ' v X c {CH^ V ) , (11) 

AH C = H C (D, H x ) = -B H D-^ V Y C {CH X D^I V ). (12) 

The scaling functions X c , Y c are universal, while B?, Bjj and C are non-universal critical am- 
plitudes, which are normalized such that X c (Cw) = l + (Cw) 2 + . . . , Y c {Cw) ~ Cw+o(Cw) 3 . 
Note that both functions are analytic for w — > 0, and AT C should be an even function of 
H x and hence w, while AH C must be an odd function of H x . From these considerations, for 
small Hi Eqs. (jl]), (§) result, rememberingil that fi c (D) — /U c (oo) = —2H C . 

An alternative argument for Eq. (|l|), which also elucidates how this equation fits together 
with the Kelvin equation (H c oc —H\/D) in the critical region for large enough D, derives 
from a consideration of phase coexistence for temperatures slightly below T C (D). If H± = 
H — 0, we would have two coexisting phases with magnetization profiles m + (z) and m~(z) = 
—m + (z) across the film, and both states have the same free energy F + (0, 0) = -F_(0, 0). Since 
these profiles are smooth functions of H and Hi, an expansion of the free energies around 
F + (0,0) {or F_(0,0), respectively} yields 

- AF+ = F+(0, 0) - F+(H, Hi) = 7^+HDL 2 + 2m+H 1 L 2 , (13) 

- AF_ = F_(0, 0) - F-(H, Hi) = ^FHDL 2 + 2rriiHiL 2 , (14) 

where m + ,m~ refer to the average over the magnetization profile in the respective states, 
and mf, rrii the layer magnetizations in the surface layer. To leading order for small Hi and 
small H in Eqs. (|T^,|T^), m + , m~ and mf, vtl{ are to be taken at zero fields, and thus satisfy 
the symmetry mr = —m + and m{ = —mf. Phase coexistence occurs for AF + = AF~ , 
and hence the Kelvin equation results, 

Hmm)x j^^n. (15) 
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Assuming then that D>(, the correlation length in the bulk, m + (z) will approach the bulk 
spontaneous magnetization m& = B^—t) 13 almost everywhere, and hence fTjZjT ps Bb(-t) 13 , 
Bfy being the respective critical amplitude. Likewise mf(D,T) approaches the surface layer 
magnetization of a semi-infinite systeml^ES vrti = B^—t)^ 1 , with B\ the corresponding 
amplitude. Therefore Eq. (|TJ) becomes in this limit 

^p.T^O^-^^Hf^, x = D\tr^oo. (16) 

u B b 

SinceS (3 « 0.325 andS-EHl fa « 0.78 - 0.80, we see that the coefficient of the term HJD 
in Eq. fll6"|) gets smaller and smaller as |t| gets smaller. Due to this vanishing coefficient in 
the limit \t\ ^0a smooth crossover between the Kelvin equation, Eqs. fllSl, [16]), and Eq. ([[]) 



becomes possible. Remembering that for D finite there is a shift of T c as given by Eq. ([|), 
we can further conclude that the critical field should be of the order 

H c (Dx, Hi) = H cocx (D, T C (D), Hi) oc - 2 Jh D -^-^ = _ 2HlD -^-^)/- ? (17) 

where in the last step the standard scaling relations (3 = 2 — a — A, flx = 2 — a — v — Ax 
were used. Eq. ( |TTD obviously is nicely compatible with Eq. (|T]). 

Eq. ( PS|) allows an interesting conclusion to be drawn on the slope of the coexistence 
curve at the critical point of the film (cf. Fig. p. We find first for the angle 8(t) describing 
this slope for T < T C (D) 

tan(0) = (3H/dT) Hl = 2 f 1 Pl " ^ {-tf^-\x = D\t\" - oo . (18) 

i c (ooJ D B b 

Since (3\ < 1 + (3, the exponent of (—t) is negative, and thus for the considered limit the 
slope diverges (i.e., varying t at very large but fixed D). However, this limit does not include 
the limiting slope at T — > T C (D) itself, since then Eqs. (|2]), ( |TT| ) yield t oc D _1 ^ u , and hence 

tan(0) oc HtD-fr-W"- 1 = HiD-^-^'". (19) 

In Landau theory, A = 3/2, Ai = 1/2 and hence the power of D vanishes, i.e., the slope is 
nonzero and finite at T C (D) in the limit D — ► oo. For the three-dimensional Ising model, the 
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best available exponent estimateSI^EHI imply (A — Ai - \)/v « 0.12 - 0.16, i.e. 9 -> 
for .D — > oo! This result also implies that for capillary condensation field mixing effectsU 
are asymptotically not very important. 

We now briefly consider the extension of the scaling theory to include finite-size effects 
due to the finite lateral linear dimension L {in Eqs. ([T^, [Hp we have assumed the limit 
L — > oo throughout}. This can simply be done by including the aspect ratio L/D as an 
additional scaling variable in Eqs. (|| which we then rewrite as follows 

U ng (D,T,H,Hi,L) « D~ 3 f(D 1 ^ v t, L/D, HD A ^ V , H x D Al ^ v ). (20) 

Since for finite L the free energy and its derivatives are smooth functions of t, it is more 
convenient to use D x l v t rather than x = D\t\ u as a scaling variable. From Eq. (|20"D , we 
immediately obtain the following scaling results for the specific heat, the magnetization and 



the susceptibility of the thin film 

C v = D a/v C{D l/v t, L/D, HD A/u , H^ 1 ^), (21) 

m = D- plv m(D l/v t, L/D, HD A/u , H^ 1 ^), (22) 

x = D~< /u x(D 1/u t, L/D, HD A/u , H^ 1 '"), (23) 



where C, m and % are appropriate scaling functions. Since we choose LL\ fixed, D fixed, 
HiD Al/ " 1, and H is chosen according to Eqs. fllB] , |I?D {in practice this is done by 
applying the reweighting technique and the equal area rule, cf. Fig. [3]} the last two arguments 
HD A I V , H X D A ^ V in Eqs. ©-(gD can be ignored in the following discussion. 

Now in the limit L — ► oo we expect that C v exhibits a logarithmic singularity for T — > 
T C (D), while the critical part of the magnetization m crit = m — m(T c (D) , H , Hi) should 
behave aS 

m crit oc [1 - T/T c (D)f\ (3 2 = 1/8, (24) 

and the susceptibility 
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X cx|l-T/T c OD)r 2 , 72 = 7/4. (25) 

For finite L, however, these singularities are all rounded off and we rather expect that 

both C v and x exhibit maxima of finite height at temperatures T^ ax (£)), T* ax (D). From 

Eqs. (gj), we readily predict (in the limit {H^D^ < 1) 

TLAD)-T c (D) Dl/u = 
i c (oo) 

r c (ooj 

with AT^ ax (L/D), AT* ax (L/D) being suitable scaling functions that describe the shift of 
these maxima positions as functions of the aspect ratio L/D. From this analysis one also 
can predictH how the height of the maxima should depend on D and L, for L >> D 

C7 ax <xD a/u ln(L/D), (28) 

X max oc D^ /v - 7/A L 7/ \ (29) 

and how the absolute value of the order paramter should decrease at T C (D), 

(\m clii \) T=Tc(D) oc D^h L -V8 . (30) 

Finally, in the limit L — > oo the D-dependence of the critical amplitudes associated with 
the two-dimensional critical behavior (see Ref. [58] for a more detailed discussion), defining 
now t = [T — T c (£>)] /T c (oo), can be read off from the following equations, 

C v ocD a/u \n\t\ , (31) 
m crit oc D^-W (-t) 1/8 , (32) 

x 0,^(7-7/4)/,^! -7/4 (33) 

Due to the crossover scaling between two- and three-dimensional critical behavior, that 
Eqs. Q20D, (|2~1"D, (]2"2"|), and describe, a singular dependence of the various critical ampli- 



tudes on film thickness results at the capillary condensation critical point. 
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IV. NUMERICAL RESULTS ON T c (D) AND H c (D) 



For locating critical points in the bulk, a convenient method is to study the fourth order 
cumulant of the order parameter 



U L (T) = 1 — (mU) I [3 (m 2 crit ) j (34) 

for a range of linear dimensions L as a function of temperature, and to look for a common 
intersection pointS {which for the universality class of the two-dimensional Ising model, 
should have the valued U* = Ul{T = T c ) = 0.610690 (1)}. In our case, we have to follow 
a path along H = H coex (T) {as shown in Fig. H(b)} when we record these cumulants, and 
since this path is not exactly known but only within some numerical error, it is clear that 
this method is more difficult to apply than for ordinary bulk Ising models. In addition, even 
for small D the data are plagued by crossover scaling effects (Fig. the curves for the 
values of L that are practically available do not intersect in a common point, but rather 
the intersection points are scattered and fall below the theoretical value U*. This failure of 
verifying the common intersection points is not unexpected, since Fig. |] includes data for 
which the aspect ratio L/D is as small as 4 (a) or even 2 (b), rather than only data for which 
L/D ^> 1. In fact, from the treatment of the previous section we can readily conclude that 

U L (T = T C (D)) = U(L/D) (35) 

and only in the limit L/D — > 00 shall we have U (00) = U*. 

An alternative and widely used recipe to find the critical temperature is to try an extrap- 
olation of the maxima of the specific heat and susceptibility versus L _1 /^or of the cumulant 
intersection points. Considering the intersection of Ul (T) and Ubi (t) with a scale factor 
b > 1, it can be argued!! that corrections to finite-size scaling lead to a shift of the inter- 
section point that varies with b proportional to b x l v — 1 for large b. Fig. [| shows some 
attempts to carry out such extrapolations, again for D = 8 and D = 28 (data for all other 
choices of D can be found inEl and look similar). These figures show that T^ ax approaches 
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T c (D) in a non-monotonic fashion, and also the curve T^ ax vs. L~ l l u is distinctly non-linear. 
Fitting asymptotic straight lines to both data sets one obtains results for T c (D) that are 
roughly compatible with each other, and with the (linear) extrapolation of the cumulant 
intersections. Although the accuracy of T c (D) obtained in this way is several orders of 
magnitude less than in the case of the bulk three-dimensional Ising modelEl, the data are 
accurate enough to allow a useful test of Eqs. ([[]), (Q). 

The consistency of our analysis can be checked further by testing for the scaling behavior 
predicted in Eqs. (|26|), (|27|), see Fig. ||. Here all data points are included for all values of 
D and L that have been studied and T C (D) is chosen such that the best data collapse is 
achieved. It is seen that the non-monotonic variation of the temperature at which the specific 
heat has its maximum is an intrinsic property of this scaling function describing the system 
shape effects in terms of the aspect ratio D/L of the simulation box. The interpolating 
curves are simple parabolic fits which translate back into the solid lines in the left part of 
Fig. H 

The values of T c (D) that we have determined as shown in Figs. [|, ^| are collected in 
Table I, which includes also our estimates for H c (D). Log- log plots of these data versus D 
almost look like straight lines, however, there is a slight but systematic curvature, and if this 
curvature were disregarded and straight lines were fitted to all the data nevertheless, the 
resulting effective exponents would systematically deviate from the theoretial predictions in 
Eqs. flf, (@). 

Better results are obtained if one fits effective exponents from successive thicknesses 
only (.0 = 4,8, 12; D = 8, 12, 16; ...;£> = 24, 28, 32), which can be extrapolated vs. 1/D 
reasonably well, and converge nicely towards the theoretical predictions (Fig. [5]), namely 
—1.587 and — (A — Ai) jv « —1.75. Conversely, if Fig. 0(b) was taken as an inde- 
pendent estimation of the exponent Ai, we would obtain Ai = 0.459(13), which indeed is 
compatible with the existing recent estimates!^ 
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V. A TEST OF TWO-DIMENSIONAL CRITICAL BEHAVIOR 



In this section, we are concerned with the question to what extent the data provide some 
evidence for the prediction (Eq. |^) that the capillary condensation critical point displays 
critical exponents of the two-dimensional Ising universality class. Since the accessible values 
of the lateral linear dimension L are not very large, however, we cannot expect that a 
regime can be reached where the parallel correlation length £| satisfies the criterion D <C 
£|| < L - only in such a regime a direct observation of these power laws would be possible. 
Hence we attempt to study the critical behavior again via a finite-size scaling analysis, 
using&B@Bi§§, 



(|m crit |) L° = M (L u t) , (36) 

X L- w = x{L u t) , (37) 
where the exponents u,v,w should take the values 

u = l/u 2 = 1, v = p/v* = 1/8, w = l2 jv 2 = 7/4. (38) 



Eqs. (|3T)|), Q3"7|), and ( |3"5D are appropriate if D -C £\\ still holds but £|| and L are of the same 
order. In practice, however, also the condition D <C £|| is hard to satisfy since we wish to 
include some data for which L/D is not very large. It then helps to relax the theoretical 
condition, Eq. (|38|), and rather treat u, v, w as effective exponents^: in this way, one can take 
into account to some extent the corrections to finite-size scaling arising from the crossover 
between two- and three-dimensional Ising critical behavior. 

Fig. |8] shows that this procedure works reasonably well, and Table II gives a listing of 
the fit exponents u,v,w, and corresponding effective exponents z/ eff — 1/u, /3 e s = v/u and 
7efr — w/u. It is seen from Table II that both u,v and w gradually change from the two- 
dimensional values towards the three-dimensional ones, although even for D = 32 one is still 
far away from the theoretical values for the latter. While /3 e g has increased significantly, 7 e fr 
within the accuracy of this estimation has hardly changed at all. If we consider an effective 
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dimensionality from the hyperscaling relation!!, defined as d e & = (7 e g + 2/5 e ff) / Teff = w + 2v, 
we find c? e g = 2.0±0.15, and there is no systematic trend with D. While the latter observation 
is in accord with a previous study using "neutral walls''^, where H C (D) = 0, we have 
obtained in the present work a much better evidence that for small D the behavior is 
compatible with two-dimensional Ising criticality than was possible in the latter modelS. 
Note also that in the present study there is a rather broad range of D where u e g > 1, which 
was not the case inS Due to the systematic problems of fitting several effective exponents 
from somewhat noisy data and the restricted range over which this fit is applicable we do 
not think that these discrepancies are a proof of non-universal crossover behavior, however. 
We feel that this problem needs a more careful study. 



VI. CONCLUSIONS 

In this paper Monte Carlo simulations have been presented attempting to test theoret- 
ical predictions about the capillary condensation critical point. Using an extension of the 
Swendsen-Wang cluster algorithm including competing surface and bulk magnetic fields, for 
Ising films with thicknesses D = 4,8, 12, 16, 24, 28 and 32 the critical temperature T C (D) 
and corresponding critical field H C (D) for a surface magnetic field Hi have been estimated. 
The data are compatible with the power laws presented about 20 years ago by Fisher and 
Nakanishi. Also the expected two-dimensional critical behavior is compatible with our data, 
though the accuracy of the resulting effective exponents is rather low (Table II) and hence 
a more convincing proof would be desirable, but is not feasible with the present computer 
ressources. 

A challenging problem that we have not solved is the development of an efficient ver- 
sion of the cluster algorithm that allows to work with surface and bulk fields that are 
not extremely small. The algorithm that we have used was much less efficient even for 
Hi = — 0.015J than for Hi = 0, and a study of capillary condensation critical points over 
the range where (H 1 / J)D Al ^ u is not small, and hence the nonlinear part of the scaling func- 
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tion Y c (CH 1 D Al / u ) would be probed, turned out not to be feasible either. Thus, in spite of 
a longstanding effort to deal with theory and simulation of capillary condensation there re- 
main still some missing links. A particularly intriguing problem is to elucidate the crossover 
between three-dimensional and two-dimensional critical behavior in these thin films. Fi- 
nally, it is hoped that the present study provides an incentive to address this problem also 
by suitable experiments. 
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TABLE I: Critical temperatures and fields 
D T C (D) H C (D)/J 



4 3.8705(3) 0.006644(32) 
8 4.2409(3) 0.002528(14) 
12 4.3561(3) 0.001367(10) 
16 4.4084(3) 0.000867(6) 
24 4.4549(5) 0.000448(3) 
28 4.4665(4) 0.000348(3) 
32 4.4749(5) 0.000279(3) 
oo 4.51152(2) 

TABLE II: Effective exponents for order parameter and susceptibility 

D U V W P eS 7eff V eS 

2- dim 1 1/8 7/4 1/8 7/4 1 

4 0.956 0.126 1.67 0.132 1.75 1.064 

8 1.018 0.136 1.72 0.133 1.69 0.982 

12 0.944 0.138 1.67 0.146 1.77 1.059 

16 0.938 0.139 1.61 0.148 1.72 1.066 

24 0.898 0.145 1.53 0.161 1.70 1.114 

28 0.853 0.141 1.48 0.165 1.74 1.172 

32 0.884 0.155 1.54 0.175 1.74 1.131 

3- dim 1.587 0.518 1.96 0.327 1.24 0.630 
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FIGURES 

FIG. 1. Schematic phase boundary for an Ising film of thickness D, where on both surfaces a 
field Hi acts, in the plane of variables temperature T and bulk field H. 

FIG. 2. Error Am as calculated from Eq. (^j plotted vs. block length iVj, for the case D = 32, 
L = 128, and two choices of Hi, Hi = (a) and H\ = —0.015 J (b). In both cases the chosen 
temperature (and bulk field H in the case of (b)) are adjusted such that the system is precisely 
at the critical point. Upper curve in each panel represents the Metropolis algorithm, lower curve 
represents the cluster algorithm. 

FIG. 3. (a) Unnormalized histogram P(m,mi) of the system with D = 32, L = 96, 
Hi/ J = —0.015, H/J = 0.00028 at k B T/J = 4.471, which is a state close to the two-phase 
coexistence line, (b) Two-phase coexistence line in the plane of variables H/J and k-oT/J for 
D = 28, estimated separately for four different choices of L from the "equal weight" -rule, showing 
also the statistical errors as estimated from Jackknife procedures!!. The two vertical lines show 
the error interval of the critical temperature. 

FIG. 4. Cumulants Ul(T) plotted vs. T for D = 8 (a) and D = 28 (b), for various choices of L 
as indicated in the figures. Dotted horizontal straight lines indicate the theoretical value U* taken 
from@. 

FIG. 5. Temperatures of specific heat and susceptibility maxima (left part) plotted vs. L~ l l v , 
and temperatures of cumulant intersections plotted vs. {b 1 ^ — 1) _1 (right part), for D = 8 (a) 
and D = 28 (b). In the left part the dashed curves show straight line fits and the solid curves 
correspond to the master curves in Fig. ^. 

FIG. 6. Master curves for temperature of the susceptibility maxima (upper part) and specific 
heat maxima (lower part) plotted vs. the inverse aspect ratio. 

FIG. 7. Plot of — 1/fcff (a) and — [(A — Ai)/v] eS (b) vs. 1/D (effective exponents were fitted 
from three successive values of D, cf. text). 
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FIG. 8. (a) Finite-size scaling plot for the critical part of the magnetization, (\m cr it\), for 
D = 4 and four choices of L as indicated, using T C (D) as quoted in Table I, and effective exponents 
u = 0.956, v = 0.126. The straight line has a slope indicating the exponent = 1/8. (b) Same 
as (a) but for D = 32, using now u = 0.884, v = 0.155. (c) Same as (b) but for the susceptibility, 
using u = 0.884, w = 1.55. The straight line indicates the exponent 72 = 7/4. 
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